7-7 mMRC breathlessness scale

Analyses of the odified Medical Research Council breathlessness scale at 28 days.

Authors
Affiliations

James Totterdell

University of Sydney

Rob Mahar

University of Melbourne

Published

July 26, 2023

Load packages 
library(ASCOTr)
library(tidyverse)
library(lubridate)
library(kableExtra)
library(patchwork)
library(cmdstanr)
library(posterior)
library(bayesplot)
library(ggdist)
library(lme4)
library(broom)
library(broom.mixed)
library(bayestestR)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))

bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))

color_scheme_set("red")
options(digits = 4)
Prepare datasets 
all_data <- readRDS(file.path(ASCOT_DATA, "all_data.rds"))

breath_labels <-c("Only breathless with strenuous exercise",
                  "Short of breath when hurrying on level ground or walking up a slight hill",
                  "Walks slower than most people of the same age because of breathlessness on level",
                  "Stops for breath after walking about 100 metres or after a few minutes on level ground",
                  "Too breathless to leave the house, or breathless when dressing or undressing")

breath_labels_short <-c("With exsercise",
                     "Up a slight hill",
                     "Slow for age",
                     "After 100 metres",
                     "Can't leave house")

# FAS-ITT
fas_itt_dat <- ASCOTr:::make_fas_itt_set(all_data) %>%
  mutate(out_mmrc_lab = fct_explicit_na(
      factor(out_mmrc_scale, 
             levels = c(1:5, 0),
             labels = c(breath_labels_short, "Not asked"))))    
fas_itt_nona_dat <- fas_itt_dat %>%
  filter(!is.na(out_mmrc_scale))

# ACS-ITT
acs_itt_dat <- ASCOTr:::make_acs_itt_set(all_data) %>%
  mutate(out_mmrc_lab = fct_explicit_na(
      factor(out_mmrc_scale, 
             levels = c(1:5, 0),
             labels = c(breath_labels_short, "Not asked"))))    
acs_itt_nona_dat <- acs_itt_dat %>%
  filter(!is.na(out_mmrc_scale))

# AVS-ITT
avs_itt_dat <- ASCOTr:::make_avs_itt_set(all_data) %>%
  mutate(out_mmrc_lab = fct_explicit_na(
      factor(out_mmrc_scale, 
             levels = c(1:5, 0),
             labels = c(breath_labels_short, "Not asked"))))    
avs_itt_nona_dat <- avs_itt_dat %>%
  filter(!is.na(out_mmrc_scale))
Load models 
ordmod0 <- compile_cmdstanr_mod(
  file.path("ordinal", "logistic_cumulative"), dir = "stan")
ordmod <- compile_cmdstanr_mod(
  file.path("ordinal", "logistic_cumulative_epoch_site"), dir = "stan")
logistic <- compile_cmdstanr_mod(
  file.path("binary", "logistic_site_epoch"), dir = "stan")
Helper functions 
make_summary_table_anticoagulation <- function(dat, format = "html") {
  tdat <- dat %>%
  group_by(
    CAssignment = factor(CAssignment, levels = c("C0", "C1", "C2", "C3", "C4"), labels = intervention_labels2()$CAssignment)) %>%
  summarise(
    Patients = n(),
    Known = sum(!is.na(out_mmrc_scale)),
        `With exercise` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 1, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 1, na.rm = TRUE)),
        `Up a slight hill` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 2, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 2, na.rm = TRUE)),
        `Slow for age` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 3, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 3, na.rm = TRUE)),
        `After 100 metres` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 4, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 4, na.rm = TRUE)),
        `Can't leave house` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 5, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 5, na.rm = TRUE)),
    `mMRC scale, Median (Q1, Q3)` = sprintf(
      "%.0f (%.0f, %.0f)", 
      median(out_mmrc_scale, na.rm = T), 
      quantile(out_mmrc_scale, 0.25, na.rm = TRUE), 
      quantile(out_mmrc_scale, 0.75, na.rm = TRUE))
  ) %>%
  bind_rows(
    dat %>%
  group_by(CAssignment = "Overall") %>%
  summarise(
    Patients = n(),
    Known = sum(!is.na(out_mmrc_scale)),
        `With exercise` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 1, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 1, na.rm = TRUE)),
        `Up a slight hill` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 2, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 2, na.rm = TRUE)),
        `Slow for age` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 3, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 3, na.rm = TRUE)),
        `After 100 metres` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 4, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 4, na.rm = TRUE)),
        `Can't leave house` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 5, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 5, na.rm = TRUE)),
    `mMRC scale, Median (Q1, Q3)` = sprintf(
      "%.0f (%.0f, %.0f)", 
      median(out_mmrc_scale, na.rm = T), 
      quantile(out_mmrc_scale, 0.25, na.rm = TRUE), 
      quantile(out_mmrc_scale, 0.75, na.rm = TRUE))
  )
  ) %>%
  rename(`Anticoagulation\nintervention` = CAssignment)
  kable(
    tdat,
    format = format,
    align = "lrrrrr",
    booktabs = TRUE,
    linesep = ""
  ) %>%
    kable_styling(
      font_size = 9,
      latex_options = "HOLD_position"
    ) %>%
    row_spec(nrow(tdat), bold = T)
}

make_summary_table_antiviral <- function(dat, format = "html") {
  tdat <- dat %>%
  group_by(
    AAssignment = factor(AAssignment, levels = c("A0", "A1", "A2"), labels = intervention_labels2()$AAssignment)) %>%
  summarise(
    Patients = n(),
    Known = sum(!is.na(out_mmrc_scale)),
        `With exercise` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 1, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 1, na.rm = TRUE)),
        `Up a slight hill` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 2, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 2, na.rm = TRUE)),
        `Slow for age` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 3, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 3, na.rm = TRUE)),
        `After 100 metres` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 4, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 4, na.rm = TRUE)),
        `Can't leave house` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 5, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 5, na.rm = TRUE)),
    `mMRC scale, Median (Q1, Q3)` = sprintf(
      "%.0f (%.0f, %.0f)", 
      median(out_mmrc_scale, na.rm = T), 
      quantile(out_mmrc_scale, 0.25, na.rm = TRUE), 
      quantile(out_mmrc_scale, 0.75, na.rm = TRUE))
  ) %>%
  bind_rows(
    dat %>%
  group_by(AAssignment = "Overall") %>%
  summarise(
    Patients = n(),
    Known = sum(!is.na(out_mmrc_scale)),
        `With exercise` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 1, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 1, na.rm = TRUE)),
        `Up a slight hill` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 2, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 2, na.rm = TRUE)),
        `Slow for age` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 3, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 3, na.rm = TRUE)),
        `After 100 metres` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 4, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 4, na.rm = TRUE)),
        `Can't leave house` = sprintf("%i (%.0f%%)", 
                                  sum(out_mmrc_scale == 5, na.rm = TRUE), 
                           100 * mean(out_mmrc_scale == 5, na.rm = TRUE)),
    `mMRC scale, Median (Q1, Q3)` = sprintf(
      "%.0f (%.0f, %.0f)", 
      median(out_mmrc_scale, na.rm = T), 
      quantile(out_mmrc_scale, 0.25, na.rm = TRUE), 
      quantile(out_mmrc_scale, 0.75, na.rm = TRUE))
  )
  ) %>%
  rename(`Antiviral\nintervention` = AAssignment)
  kable(
    tdat,
    format = format,
    align = "lrrrrr",
    booktabs = TRUE,
    linesep = ""
  ) %>%
    kable_styling(
      font_size = 9,
      latex_options = "HOLD_position"
    ) %>%
    row_spec(nrow(tdat), bold = T)
}

Descriptive

Anticoagulation

Summary of mMRC outcome by arm 
save_tex_table(
  make_summary_table_anticoagulation(acs_itt_dat, "latex"), 
  file.path("outcomes", "secondary", "7-7-anticoagulation-summary"))
make_summary_table_anticoagulation(acs_itt_dat)
Anticoagulation intervention Patients Known With exercise Up a slight hill Slow for age After 100 metres Can't leave house mMRC scale, Median (Q1, Q3)
Low-dose 610 577 40 (7%) 47 (8%) 12 (2%) 14 (2%) 2 (0%) 0 (0, 0)
Intermediate-dose 613 584 50 (9%) 39 (7%) 11 (2%) 8 (1%) 2 (0%) 0 (0, 0)
Low-dose with aspirin 283 271 26 (10%) 16 (6%) 15 (6%) 1 (0%) 1 (0%) 0 (0, 0)
Therapeutic-dose 50 44 2 (5%) 4 (9%) 3 (7%) 2 (5%) 0 (0%) 0 (0, 0)
Overall 1556 1476 118 (8%) 106 (7%) 41 (3%) 25 (2%) 5 (0%) 0 (0, 0)
Table 1: Summary of mMRC scale at day 28 by treatment group, anticoagulation domain, FAS-ITT.
Code 
p1 <- fas_itt_nona_dat %>%
  group_by(CAssignment = factor(CAssignment, labels = str_replace(intervention_labels()$CAssignment, "<br>", "\n"))) %>%
  summarise(n = n()) %>%
  ggplot(., aes(CAssignment, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Anticoagulation intervention")
p2 <- fas_itt_nona_dat %>%
  filter(!is.na(out_mmrc_scale)) %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      labels = str_replace(intervention_labels()$CAssignment, "<br>", "\n")), 
    out_mmrc_scale = ordered(as.integer(out_mmrc_scale), levels = 0:5)
  ) %>% 
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  ggplot(., aes(CAssignment, p)) +
  geom_col(aes(fill = out_mmrc_scale)) +
  scale_fill_viridis_d(option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  labs(fill = "mMRCbs", y = "Cumulative proportion", x = "Anticoagulation intervention")
p <- p1 | p2
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "7-7-anticoagulation.pdf"), p2, height = 3, width = 6)
p2

Antiviral

Summary of mMRC outcome by arm 
save_tex_table(
  make_summary_table_antiviral(avs_itt_dat, "latex"), 
  file.path("outcomes", "secondary", "7-7-antiviral-summary"))
make_summary_table_antiviral(avs_itt_dat)
Antiviral intervention Patients Known With exercise Up a slight hill Slow for age After 100 metres Can't leave house mMRC scale, Median (Q1, Q3)
Standard of care 73 67 5 (7%) 15 (22%) 5 (7%) 9 (13%) 1 (1%) 1 (0, 2)
Nafamostat 82 73 9 (12%) 12 (16%) 5 (7%) 8 (11%) 2 (3%) 0 (0, 2)
Overall 155 140 14 (10%) 27 (19%) 10 (7%) 17 (12%) 3 (2%) 1 (0, 2)
Table 2: Summary of mMRC scale at day 28 by treatment group, antiviral domain, FAS-ITT.
Code 
p1 <- fas_itt_nona_dat %>%
  group_by(AAssignment = factor(AAssignment, labels = str_replace(intervention_labels()$AAssignment, "<br>", "\n"))) %>%
  summarise(n = n()) %>%
  ggplot(., aes(AAssignment, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Antiviral intervention")
p2 <- fas_itt_nona_dat %>%
  filter(!is.na(out_mmrc_scale)) %>%
  dplyr::count(
    AAssignment = factor(
      AAssignment, 
      labels = str_replace(intervention_labels()$AAssignment, "<br>", "\n")), 
    out_mmrc_scale = ordered(as.integer(out_mmrc_scale), levels = 0:5)
  ) %>%
  group_by(AAssignment) %>%
  mutate(p = n / sum(n)) %>%
  ggplot(., aes(AAssignment, p)) +
  geom_col(aes(fill = out_mmrc_scale)) +
  scale_fill_viridis_d(option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  labs(fill = "mMRCbs", y = "Cumulative proportion", x = "Antiviral intervention")
p <- p1 | p2
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "7-7-antiviral.pdf"), p2, height = 3, width = 6)
p

Age

Code 
pdat <- fas_itt_nona_dat %>%
  dplyr::count(
    agegrp = cut(AgeAtEntry, c(18, seq(25, 75, 5), 100), include.lowest = T, right = F),
    mrc = ordered(out_mmrc_scale, levels = 0:5), 
    .drop = F) %>%
  group_by(agegrp) %>%
  mutate(p = n / sum(n))
pdat2 <- pdat %>%
  group_by(agegrp) %>%
  summarise(n = sum(n))
p1 <- ggplot(pdat2, aes(agegrp, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry") +
  geom_vline(xintercept = 8.5, linetype = 2) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45))
p2 <- ggplot(pdat, aes(agegrp, p)) +
  geom_col(aes(fill = mrc)) +
  labs(x = "Age", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("mMRCbs", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  geom_vline(xintercept = 8.5, linetype = 2) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
  theme(legend.key.size = unit(0.5, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-7-age.pdf"), p, height = 2.5, width = 6)
p

Figure 1: Distribution of mMRC scale at day 28 by age group, FAS-ITT

Sex

Code 
pdat <- fas_itt_nona_dat %>%
  dplyr::count(
    Sex,
    mrc = ordered(out_mmrc_scale, levels = 0:5), 
    .drop = F) %>%
  group_by(Sex) %>%
  mutate(p = n / sum(n))
pdat2 <- pdat %>%
  group_by(Sex) %>%
  summarise(n = sum(n))
p1 <- ggplot(pdat2, aes(Sex, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  labs(y = "Number of\nparticipants",
       x = "Sex")
p2 <- ggplot(pdat, aes(Sex, p)) +
  geom_col(aes(fill = mrc)) +
  labs(x = "Sex", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("mMRCbs", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  theme(legend.key.size = unit(0.5, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-7-sex.pdf"), p, height = 2.5, width = 6)
p

Figure 2: Distribution of mMRC scale at day 28 by sex, FAS-ITT

Oxygen

Code 
pdat <- fas_itt_nona_dat %>%
  dplyr::count(
    supp_oxy2 = fct_explicit_na(factor(supp_oxy2, levels = 0:1, labels = c("Not required", "Required"))),
    mrc = ordered(out_mmrc_scale, levels = 0:5), 
    .drop = F) %>%
  group_by(supp_oxy2) %>%
  mutate(p = n / sum(n))
pdat2 <- pdat %>%
  group_by(supp_oxy2) %>%
  summarise(n = sum(n))
p1 <- ggplot(pdat2, aes(supp_oxy2, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  labs(y = "Number of\nparticipants",
       x = "Supplemental oxygen")
p2 <- ggplot(pdat, aes(supp_oxy2, p)) +
  geom_col(aes(fill = mrc)) +
  labs(x = "Supplemental oxygen", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("mMRCbs", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  theme(legend.key.size = unit(0.5, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-7-oxygen.pdf"), p, height = 2.5, width = 6)
p

Figure 3: Distribution of mMRC scale at day 28 by supplemental oxygen requirement, FAS-ITT

Country

Code 
pdat <- fas_itt_nona_dat %>%
  dplyr::count(
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand"),
                     labels = c("India", "Australia", "Nepal", "New\nZealand")),
    mrc = ordered(out_mmrc_scale, levels = 0:5), 
    .drop = F) %>%
  group_by(Country) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
pdat2 <- pdat %>%
  group_by(Country) %>%
  summarise(n = sum(n))
p1 <- ggplot(pdat2, aes(Country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")
p2 <- ggplot(pdat, aes(Country, p)) +
  geom_col(aes(fill = mrc)) +
  labs(x = "Country", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("mMRCbs", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  theme(legend.key.size = unit(0.5, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-7-country.pdf"), p, height = 2.5, width = 6)
p

Figure 4: Distribution of mMRC scale at day 28 by country, FAS-ITT.

Site

Code 
pdat <- fas_itt_nona_dat %>%
  dplyr::count(
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand"),
                     labels = c("India", "Australia", "Nepal", "New\nZealand")),
    Site = fct_infreq(Location),
    mrc = ordered(out_mmrc_scale, levels = 0:5)) %>%
  complete(mrc = ordered(0:5), nesting(Country, Site), fill = list(n = 0)) %>%
  group_by(Country, Site) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup() %>%
  mutate(
    Country = droplevels(Country),
    Site = droplevels(Site)
  )
pdat2 <- pdat %>%
  group_by(Country, Site) %>%
  summarise(n = sum(n)) %>%
  ungroup()
p1 <- ggplot(pdat2, aes(Site, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(pdat, aes(Site, p)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col(aes(fill = mrc)) +
  labs(x = "Site", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("mMRCbs", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F, ncol = 1)) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
  theme(legend.key.size = unit(0.25, "lines"))
p <- p1 / p2 +
  plot_layout(guides = 'collect')
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-7-country-site.pdf"), p, height = 4, width = 6.25)
p

Figure 5: Distribution of mMRC scale by study site, FAS-ITT.

Calendar Time

Code 
pdat <- fas_itt_nona_dat %>%
  filter(!is.na(out_mmrc_scale)) %>%
  dplyr::count(
    yr = year(RandDate), mth = month(RandDate),
    mrc = ordered(out_mmrc_scale, levels = 0:5), 
    .drop = F) %>%
  group_by(yr, mth) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p1 <- pdat %>%
  group_by(yr, mth) %>%
  summarise(n = sum(n)) %>%
  ggplot(., aes(mth, n))  +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p2 <- ggplot(pdat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
  geom_col(aes(fill = mrc)) +
  labs(x = "Calendar date (month of year)", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("mMRCbs", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  theme(legend.key.size = unit(0.5, "lines")) +
  scale_x_continuous(breaks = 1:12)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-7-calendar-time.pdf"), p, height = 2, width = 6)
p

Figure 6: Distribution of mMRC scale by calendar time, FAS-ITT.

Sample Cumulative Logits

Proportional odds looks reasonable for all logits.

Plot sample cumulative logits 
trt_counts <- fas_itt_nona_dat %>%
  dplyr::count(CAssignment, out_mmrc_scale) %>%
  complete(CAssignment, out_mmrc_scale, fill = list(n = 0)) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n))
trt_logit <- trt_counts %>% 
  group_by(CAssignment) %>% 
  mutate(clogit = logit(cumsum(p))) %>%
  group_by(out_mmrc_scale) %>%
  mutate(rel_clogit = clogit - mean(clogit)) %>%
ggplot(., aes(out_mmrc_scale, rel_clogit)) +
  facet_wrap( ~ CAssignment) +
  geom_point() +
  labs(y = "Relative (to mean) sample cumulative logit")
trt_logit

Inspect sample cumulative logits, anticoagulation.

Some indication that a contrained partial proportional odds model might have better ‘fit’, but sample size is very small, so sticking with PO model for consistency among results.

Plot sample cumulative logits 
trt_counts <- fas_itt_nona_dat %>%
  dplyr::count(AAssignment, out_mmrc_scale) %>%
  complete(AAssignment, out_mmrc_scale, fill = list(n = 0)) %>%
  group_by(AAssignment) %>%
  mutate(p = n / sum(n))
trt_logit <- trt_counts %>% 
  group_by(AAssignment) %>% 
  mutate(clogit = logit(cumsum(p))) %>%
  group_by(out_mmrc_scale) %>%
  mutate(rel_clogit = clogit - mean(clogit)) %>%
ggplot(., aes(out_mmrc_scale, rel_clogit)) +
  facet_wrap( ~ AAssignment) +
  geom_point() +
  labs(y = "Relative (to mean) sample cumulative logit")
trt_logit

Inspect sample cumulative logits, antiviral.

Modelling

FAS-ITT

Fit primary model 
res <- fit_primary_model(
  fas_itt_nona_dat %>% mutate(out_mmrc_scale = out_mmrc_scale + 1),
  ordmod,
  outcome = "out_mmrc_scale",
  intercept = FALSE
)
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Ineligible aspirin", "Age \u2265 60", "Female", "Oxygen requirement", "Australia/New Zealand", "Nepal")
Odds ratio summary table 
save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-7-primary-model-fas-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Nafamostat 0.68 (0.34, 1.37) 0.73 (0.27) 0.86
Intermediate-dose 0.74 (0.53, 1.02) 0.75 (0.13) 0.97
Low-dose with aspirin 1.11 (0.73, 1.67) 1.14 (0.24) 0.31
Therapeutic-dose 0.86 (0.34, 2.01) 0.94 (0.44) 0.64
Ineligible aspirin 1.43 (0.56, 3.44) 1.58 (0.76) 0.23
Age ≥ 60 1.98 (1.44, 2.70) 2.01 (0.32) 0.00
Female 1.03 (0.77, 1.37) 1.04 (0.15) 0.43
Oxygen requirement 1.17 (0.83, 1.64) 1.19 (0.21) 0.18
Australia/New Zealand 2.54 (0.68, 9.16) 3.16 (2.39) 0.08
Nepal 0.44 (0.10, 2.38) 0.63 (0.77) 0.85
Posterior summaries for model parameters (fixed-effects), FAS-ITT.
Odds ratio posterior densities 
p <- plot_or_densities(c(res$drws$AOR, res$drws$COR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "primary", "7-7-primary-model-fas-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Figure 7: Posterior densities for odds ratio contrasts.
Odds ratio posterior summary for epoch and site 
p <- plot_epoch_site_terms(
  res$drws$gamma_epoch,
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-7-primary-model-epoch-site-terms-fas-itt.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Summary of epoch and site posterior odds ratios.
Code 
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.7977 0.8511 0.7592 0.7486 0.8220 0.7691 0.7539 0.7562
Code 
mcmc_trace(res$drws["beta"])

Code 
mcmc_trace(res$drws["gamma_site"])

Code 
mcmc_trace(res$drws["gamma_epoch"])

Posterior Predictive

Posterior predictive figure 
y_raw <- as.integer(levels(res$dat$y_raw))
y_ppc <- res$drws$y_ppc
y_ppc_raw <- rfun(\(x) y_raw[x])(y_ppc)
ppc_dat <- bind_cols(fas_itt_nona_dat, tibble(y_ppc = y_ppc_raw))
grp_ppc1 <- function(grp) {
  ppc_dat  %>%
  group_by(grp = !!grp) %>%
  summarise(
    y_1 = mean(out_mmrc_scale == 0),
    ypp_1 = rvar_mean(y_ppc == 0),
    y_1 = mean(out_mmrc_scale == 1),
    ypp_1 = rvar_mean(y_ppc == 1),
    y_2 = mean(out_mmrc_scale <= 2),
    ypp_2 = rvar_mean(y_ppc <= 2),
    y_3 = mean(out_mmrc_scale <= 3),
    ypp_3 = rvar_mean(y_ppc <= 3),
    y_4 = mean(out_mmrc_scale <= 4),
    ypp_4 = rvar_mean(y_ppc <= 4),
    y_5 = mean(out_mmrc_scale <= 5),
    ypp_5 = rvar_mean(y_ppc <= 5)
  ) %>%
  pivot_longer(y_1:ypp_5, names_to = c("response", "mrc"), names_sep = "_", values_to = "posterior") %>%
  mutate(mrc = as.numeric(mrc))
}
grp_ppc2 <- function(grp) {
  ppc_dat  %>%
  group_by(grp = !!grp) %>%
  summarise(
    y_1 = mean(out_mmrc_scale == 0),
    ypp_1 = rvar_mean(y_ppc == 1),
    y_2 = mean(out_mmrc_scale == 1),
    ypp_2 = rvar_mean(y_ppc == 2),
    y_3 = mean(out_mmrc_scale == 2),
    ypp_3 = rvar_mean(y_ppc == 3),
    y_4 = mean(out_mmrc_scale == 3),
    ypp_4 = rvar_mean(y_ppc == 4),
    y_5 = mean(out_mmrc_scale == 4),
    ypp_5 = rvar_mean(y_ppc == 5),
    y_6 = mean(out_mmrc_scale == 5),
    ypp_6 = rvar_mean(y_ppc == 6)
  ) %>%
  pivot_longer(y_1:ypp_6, names_to = c("response", "mrc"), 
               names_sep = "_", values_to = "posterior") %>%
  mutate(mrc = as.numeric(mrc))
}

plot_grp_ppc <- function(dat) {
  ggplot(dat %>% filter(response == "ypp"), aes(x = mrc)) +
    facet_wrap( ~ grp, nrow = 1) +
    stat_slabinterval(aes(ydist = posterior))  +
    geom_point(data = dat %>% filter(response == "y"), 
               aes(x = mrc, y = mean(posterior)),
               colour = "red",
               shape = 23) +
    labs(x = "mMRC outcome", 
         y = "Probability")
}

plot_grp_ppc <- function(dat, lab = "", xlab = "Probability") {
  ggplot(dat %>% filter(response == "ypp"), aes(y = mrc)) +
    facet_wrap( ~ grp, nrow = 1) +
    stat_interval(aes(xdist = posterior), size = 2)  +
    geom_point(data = dat %>% filter(response == "y"), 
               aes(y = mrc, x = mean(posterior)),
               colour = "red",
               shape = 23) +
    scale_x_continuous(xlab, breaks = c(0, 0.5),
                       sec.axis = sec_axis(~ . , name = lab, breaks = NULL, labels = NULL)) +
    scale_y_continuous("WHO\noutcome", breaks = c(1,3,5,7)) +
    theme(strip.text = element_text(size = rel(0.7)),
          axis.title.x = element_text(size = rel(0.7)),
          axis.text.x = element_text(size = rel(0.65)),
          axis.title.y = element_text(size = rel(0.75)),
          axis.title.x.bottom = element_blank())
}

pp_epoch <- grp_ppc2(sym("epoch")) %>% 
  mutate(grp = fct_inorder(factor(grp)))
pp_A <- grp_ppc2(sym("AAssignment"))
pp_C <- grp_ppc2(sym("CAssignment"))
pp_ctry <- grp_ppc2(sym("Country"))
pp_site <- grp_ppc2(sym("site")) %>%
  left_join(ppc_dat %>% dplyr::count(site, Country), by = c("grp" = "site"))
p0 <- plot_grp_ppc(pp_A, "Antiviral", "")
p1 <- plot_grp_ppc(pp_C, "Anticoagulation", "")
p2 <- plot_grp_ppc(pp_ctry, "Country", "") 
p3 <- plot_grp_ppc(pp_epoch, "Epoch", "")
p4 <- plot_grp_ppc(pp_site %>% filter(Country == "IN"), "Sites India", "")
p5 <- plot_grp_ppc(pp_site %>% filter(Country == "AU"), "Sites Australia", "")
p6 <- plot_grp_ppc(pp_site %>% filter(Country == "NP"), "Sites Nepal", "")
p7 <- plot_grp_ppc(pp_site %>% filter(Country == "NZ"), "Sites New Zealand", "")
p <- (p0 | p1 | p2) / p3 / p4 / p5 / (p6 | p7)+
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")
pth <- file.path("outputs", "figures", "outcomes", "secondary",
                 "7-7-primary-model-fas-itt-ppc.pdf")
ggsave(pth, p, width = 6, height = 5.75)
p

Figure 8: Posterior predictive check, FAS-ITT.

ACS-ITT

Fit primary model 
res <- fit_primary_model(
  acs_itt_nona_dat %>% mutate(out_mmrc_scale = out_mmrc_scale + 1),
  ordmod,
  outcome = "out_mmrc_scale", 
  intercept = FALSE
)
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Ineligible aspirin", "Age \u2265 60", "Female", "Oxygen requirement", "Australia/New Zealand", "Nepal")
Odds ratio summary table 
save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-7-primary-model-acs-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Nafamostat 0.78 (0.34, 1.73) 0.84 (0.36) 0.73
Intermediate-dose 0.73 (0.52, 1.01) 0.74 (0.13) 0.97
Low-dose with aspirin 1.09 (0.72, 1.65) 1.12 (0.24) 0.33
Therapeutic-dose 0.86 (0.35, 2.00) 0.94 (0.44) 0.64
Ineligible aspirin 1.30 (0.48, 3.29) 1.45 (0.73) 0.29
Age ≥ 60 1.91 (1.39, 2.63) 1.94 (0.32) 0.00
Female 1.07 (0.80, 1.43) 1.08 (0.16) 0.33
Oxygen requirement 1.20 (0.85, 1.71) 1.22 (0.22) 0.15
Australia/New Zealand 2.26 (0.57, 8.35) 2.82 (2.13) 0.12
Nepal 0.45 (0.10, 2.35) 0.65 (0.68) 0.84
Posterior summaries for model parameters (fixed-effects), ACS-ITT.
Odds ratio posterior densities 
p <- plot_or_densities(c(res$drws$AOR, res$drws$COR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "primary", "7-7-primary-model-acs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Figure 9: Posterior densities for odds ratio contrasts, ACS-ITT.
Odds ratio posterior summary for epoch and site 
p <- plot_epoch_site_terms(
  res$drws$gamma_epoch,
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-7-primary-model-epoch-site-terms-acs-itt.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Figure 10: Summary of epoch and site posterior odds ratios, ACS-ITT.
Code 
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8588 0.7435 0.8051 0.7942 0.7855 0.7693 0.7807 0.7988
Code 
mcmc_trace(res$drws["beta"])

Code 
mcmc_trace(res$drws["gamma_site"])

Code 
mcmc_trace(res$drws["gamma_epoch"])

AVS-ITT

Code 
p1 <- avs_itt_nona_dat %>%
  group_by(CAssignment = factor(CAssignment, labels = str_replace(intervention_labels()$CAssignment, "<br>", "\n"))) %>%
  summarise(n = n()) %>%
  ggplot(., aes(CAssignment, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Anticoagulation intervention")
p2 <- fas_itt_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      labels = str_replace(intervention_labels()$CAssignment, "<br>", "\n")), 
    out_mmrc_scale = ordered(as.integer(out_mmrc_scale), levels = 0:5)
  ) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  ggplot(., aes(CAssignment, p)) +
  geom_col(aes(fill = out_mmrc_scale)) +
  scale_fill_viridis_d(option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = F)) +
  labs(fill = "mMRC scale", y = "Cumulative proportion", x = "Anticoagulation intervention")
p <- p1 | p2
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "7-7-anticoagulation-avs-itt.pdf"), p2, height = 3, width = 6)
p2

Code 
res <- fit_primary_model(
  avs_itt_nona_dat %>% mutate(out_mmrc_scale = out_mmrc_scale + 1),
  ordmod0,
  outcome = "out_mmrc_scale",
  vars = c("agegte60", "sexF", "supp_oxy2", "crp_tertile"),
  beta_sd_var = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5),
  intercept = FALSE
)
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Age \u2265 60", "Female", "Oxygen requirement", "CRP (2nd tertile)", "CRP (3rd tertile)", "CRP (unknown)")
Odds ratio summary table 
save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-7-primary-model-avs-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Nafamostat 0.79 (0.40, 1.58) 0.84 (0.30) 0.75
Intermediate-dose 0.91 (0.40, 2.07) 0.99 (0.44) 0.59
Low-dose with aspirin 0.57 (0.14, 2.12) 0.71 (0.54) 0.80
Therapeutic-dose 1.56 (0.55, 4.40) 1.79 (1.04) 0.20
Age ≥ 60 2.32 (1.14, 4.87) 2.50 (0.96) 0.01
Female 1.43 (0.72, 2.81) 1.51 (0.54) 0.15
Oxygen requirement 1.67 (0.79, 3.67) 1.81 (0.75) 0.10
CRP (2nd tertile) 0.77 (0.33, 1.77) 0.84 (0.38) 0.73
CRP (3rd tertile) 0.86 (0.37, 1.97) 0.94 (0.42) 0.64
CRP (unknown) 1.04 (0.33, 3.05) 1.21 (0.72) 0.47
Posterior summaries for model parameters (fixed-effects), AVS-ITT.
Odds ratio posterior densities 
p <- plot_or_densities(c(res$drws$AOR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "primary", "7-7-primary-model-avs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Figure 11: Posterior densities for odds ratio contrasts, AVS-ITT.
Code 
mcmc_trace(res$drws["beta"])

End of script.